The present invention is directed to a system and method for reduction of x-ray scatter in imaging systems and more particularly to such a system and method which use both spatial and temporal interpolation for such reduction. The present invention is usable with several imaging techniques, including fan beam and cone beam CT.
Scatter reduction and correction are required for both medical (e.g., clinical, small animal imaging) and nonmedical imagine applications (e.g., explosive detection, nondestructive testing, industrial imagine guided manufacturing applications). Compton scatter interaction always exists in the x-ray energy range useful for imaging (10 keV to a few MeV). In some segments of the x-ray spectrum, X-ray scatter interaction is dominant. For example, within the effective energy range of medical diagnostic CT (60xcx9c80 keV), it is recognized that the Compton interaction ( greater than 80%) plays a predominant role over photoelectrical interaction ( less than 10%) and coherent interaction ( less than 5%) when x-ray photons pass through water-like soft tissues, although slightly higher percentages for both of them are observed when x-ray photons pass through compact bone. However, scatter intensity detected by a detector usually contributes to noise, not to a useful signal for imaging. Therefore, scatter reduction and correction are required to improve reconstruction accuracy of linear attenuation coefficient (LAC) distribution for both medical (clinical, small animal imaging and nonmedical imagine applications (explosive detection, nondestructive testing, industrial imagine guided manufacturing applications).
Before the mechanism in which x-ray scatter interferes with cone beam volume computed tomography (CBVCT) imaging is described, it is instructive to know how x-ray scattering deteriorates the quality of projection imaging. FIG. 1A shows an example in which a disc 102 is embedded in a uniform cylindrical object 104 whose LAC is slightly larger than that of the disc 102. If the object 104 is disposed in a cone beam C emitted by a cone-beam source 106, and the cone beam C is then received by a detector 108, a projection image is acquired as shown in FIG. 1B. Due to the LAC discrepancy, the number of x-ray photons detected in the area corresponding to the disc is N+xcex94N, but that detected in the background area is N, with no x-ray photons scattered. By limiting to an additive-noise-free situation, the ability to visualize such a structure is generally evaluated by contrast, which is defined as                     C        =                                            N              +                              Δ                ⁢                                  xe2x80x83                                ⁢                N                            -              N                        N                    =                                    Δ              ⁢                              xe2x80x83                            ⁢              N                        N                                              (        1        )            
It is known that the flux of x-ray photons observes the Poisson distribution, while the transmittance of x-ray photons through an object observes the binomial distribution. The cascading of a Poisson process and a binomial process is still a Poisson process. Hence, if the mean number of the x-ray photons transmitting the disc and its surrounding area are N+xcex94N and N respectively, their corresponding standard deviations are {square root over (N+xcex94N)} and {square root over (N)}. From the perspective of system analysis, the scattered x-ray photons behave like additive noise. In digital projection imaging where the display window of an ROI (region of interest) can be adjusted arbitrarily, the ability to visualize such a local structure is more appropriately measured by SNR (signal-to-noise ratio), which is defined as                     SNR        =                                            N              +                              Δ                ⁢                                  xe2x80x83                                ⁢                N                            -              N                                      N                                =                                                    Δ                ⁢                                  xe2x80x83                                ⁢                N                                            N                                      =                          C              ⁢                              N                                                                        (        2        )            
On the other hand, letting the mean number of the scattered x-ray photons be Ns, the scatter-to-primary ratio (SPR) is defined as                     SPR        =                              N            s                    N                                    (        3        )            
and the scatter degradation factor (SDF) is defined as                     SDF        =                              N                          N              +                              N                s                                              =                                    1                              1                +                                                      N                    s                                    /                  N                                                      =                          1                              1                +                SPR                                                                        (        4        )            
Consequently, (1) and (2) respectively degrade to                               C          s                =                                            N              +                              N                s                            +                              Δ                ⁢                                  xe2x80x83                                ⁢                N                            -                              (                                  N                  +                                      N                    s                                                  )                                                    N              +                              N                s                                              =                                                    Δ                ⁢                                  xe2x80x83                                ⁢                N                                            N                +                                  N                  s                                                      =                                          C                                  1                  +                                                            N                      s                                        /                    N                                                              =                              SDF                ·                C                                                                        (        5        )                                          SNR          s                =                                            N              +                              N                s                            +                              Δ                ⁢                                  xe2x80x83                                ⁢                N                            -                              (                                  N                  +                                      N                    s                                                  )                                                                    N                +                                  N                  s                                                              =                                                    Δ                ⁢                                  xe2x80x83                                ⁢                N                                                              N                  +                                      N                    s                                                                        =                                          SNR                                                      1                    +                                                                  N                        s                                            /                      N                                                                                  =                                                SDF                                ·                SNR                                                                        (        6        )            
in the presence of x-ray scatter. In other words, due to scattered x-ray photons, the local contrast and SNR of the structure are deteriorated by the factors of SDF and {square root over (SDF)} in a projection image, respectively.
Suppose that Ip(i,j)(ixcex5I,jxcex5J) refers to the image formed by the primary x-ray photons, and Is(i,j)(ixcex5I,jxcex5J) the image formed by the scattered photons, where I and J are the vertical and horizontal dimension of a projection image. As a result, a 2D SPR distribution is defined as                               SPR          ⁡                      (                          i              ,              j                        )                          =                                                            I                s                            ⁡                              (                                  i                  ,                  j                                )                                                                    I                p                            ⁡                              (                                  i                  ,                  j                                )                                              ⁢                      xe2x80x83                    ⁢                      (                                          i                ∈                I                            ,                              j                ∈                J                                      )                                              (        7        )            
and subsequently                               SDF          ⁡                      (                          i              ,              j                        )                          =                              1                          1              +                              SPR                ⁡                                  (                                      i                    ,                    j                                    )                                                              ⁢                      xe2x80x83                    ⁢                      (                                          i                ∈                I                            ,                              j                ∈                J                                      )                                              (        8        )            
It has been found that, given an object, the distribution of Is(i,j) is dependent on its structure, thickness and field of view (FOV) through which Ip(i,j) is formed. However, regardless of how the distribution of Ip(i,j) fluctuates, the variation of Is(i,j) is so smooth that it can be approximately treated as a 2D spatial low-pass filtering of Ip(i,J), and several low-pass filtering models associated with different filter kernels have been proposed. That means that a strong spatial correlation exists between neighboring pixels. Intuitively, scatter distribution in each projection image can be recovered from its spatial samples using either interpolation methods, such as cubic spline interpolation or bi-linear interpolation methods, or a convolution operation with a selected convolution kernel (a low pass filtering in the frequency domain).
On the other hand, the SPR(i,j) or SDF(i,j) usually fluctuates very much, especially where Ip(i,j) is relatively low. Hence, the SPR(i,j) or SDF(i,j) is usually taken into account to reflect the severity of the x-ray scatter in projection imaging.
A CBVCT (cone beam volume computed tomography) image is reconstructed from a set of consecutive 2D projection images that are sequentially acquired. Intrinsically, the x-ray scatter interferes with the X-ray transform non-linearly. Once the X-ray transform is acquired, the artifact caused by the x-ray scatter in a tomographic image is reconstruction-algorithm-dependent. Thus, the experimental investigation of the x-ray scatter artifact in tomographic imaging is preferable in practice. In conventional CT (computed tomography), due to the adoption of a slit collimator, the severity of the longitudinally scattered x-ray photons is reduced to a secondary order in comparison to that of the transversely scattered ones. To further reduce the transversely scattered x-ray, other measures, such as the bow-tie x-ray attenuator, and the post-patient out-of-slice collimator used in a 3rd-generation CT or the reference detector used in a 4th-generation CT, are usually taken into account.
Unfortunately, for CBVCT, the slit collimator has to be removed to fully use the generated cone shaped x-ray beam. Hence, much more severe x-ray scatter interference is expected in CBVCT, although the air gap technique and beam-shaping (bow-tie) attenuator are still useful for reducing scatter. Supposing Ip(x,y) refers to the projection image formed by primary x-ray photons, and Is(x,y) that formed by scattered x-ray photons, we have
It(x,y)=Ip(x,y)+Is(x,y)xe2x80x83xe2x80x83(9)
As a result, a 2D SPR distribution is defined as                               SPR          ⁡                      (                          x              ,              y                        )                          =                                            I              s                        ⁡                          (                              x                ,                y                            )                                                          I              p                        ⁡                          (                              x                ,                y                            )                                                          (        10        )            
By substituting (10) into (9), one gets
It(x,y)=Ip(x,y)[1.0+SPR(x,y)]xe2x80x83xe2x80x83(11)
Furthermore, by taking the incident x-ray intensity distribution I0(x,y) into account and applying logarithm, we get the X-ray transform data for CB reconstruction                               P          ⁡                      (                          x              ,              y                        )                          =                              ln            ⁢                                                            I                  0                                ⁡                                  (                                      x                    ,                    y                                    )                                                                              I                  i                                ⁡                                  (                                      x                    ,                    y                                    )                                                              =                                    ln              ⁢                                                                    I                    0                                    ⁡                                      (                                          x                      ,                      y                                        )                                                                                                              I                      p                                        ⁡                                          (                                              x                        ,                        y                                            )                                                        ⁡                                      [                                          1.0                      +                                              SPR                        ⁡                                                  (                                                      x                            ,                            y                                                    )                                                                                      ]                                                                        =                                                            ln                  ⁢                                                                                    I                        0                                            ⁡                                              (                                                  x                          ,                          y                                                )                                                                                                            I                        p                                            ⁡                                              (                                                  x                          ,                          y                                                )                                                                                            -                                  ln                  ⁡                                      [                                          1.0                      +                                              SPR                        ⁡                                                  (                                                      x                            ,                            y                                                    )                                                                                      ]                                                              ≡                                                                    P                    p                                    ⁡                                      (                                          x                      ,                      y                                        )                                                  +                                                      P                    s                                    ⁡                                      (                                          x                      ,                      y                                        )                                                                                                          (        12        )            
where                                           P            p                    ⁡                      (                          x              ,              y                        )                          =                  ln          ⁢                                                    I                0                            ⁡                              (                                  x                  ,                  y                                )                                                                    I                p                            ⁡                              (                                  x                  ,                  y                                )                                                                        (        13        )            
are the X-ray transform distribution corresponding to primary and scattered x-ray photons respectively. Hence, due to scatter interference, the CB reconstruction of an object becomes
{circumflex over (f)}({right arrow over (r)})+xcex94f({right arrow over (r)})xe2x80x83xe2x80x83(15)
i.e., the sum of the CB reconstruction corresponding to the projection image formed by primary x-ray photons and that formed by scattered ones.
Inferred from what is presented above, it is not difficult to obtain the following important observations concerning the reconstruction error introduced by x-ray scatter (RE) in CBVCT:
Is(x,y) is the source of Re;
RE is determined by SPR distribution xe2x88x92ln[1.0+SPR(x,y)], rather than merely Is(x,y). In other words, RE is determined by not only Is(x,y) but also Ip(x,y);
RE is always negative because xe2x88x92ln[1.0+SPR(x,y)] is less than zero. This is the rationale behind the cupping effect or shading effect in CBVCT images;
While SPR(x,y) greater than  greater than 1, i.e., the x-ray scatter interference is extraordinarily severe, RE is equivalently proportional to xe2x88x92ln SPR(x,y);
While SPR(x,y) less than  less than 1, i.e., the x-ray scatter interference is moderate, RE is equivalently proportional to xe2x88x92SPR(x,y);
The RE caused by x-ray scatter interference can be alleviated through two approaches: (a) reducing Is(x,y); (b) increasing Ip(x,y); and
It should be emphasized that a significant drop in intensity of Ip(x,y) formed by primary x-ray photons is to be avoided in practice because the resultant increase in xe2x88x92ln[1.0+SPR(x,y)] can be extremely high and results in severe RE in CBVCT.
The simplest approach to suppress x-ray scatter is the air gap, which has been employed in projection imaging, particularly chest radiography, for many years. A CBVCT inherits the air gap approach due to the distance maintained between x-ray detector and the object to be scanned. As illustrated in FIG. 2, T represents the thickness of the object 202 to be imaged, xp the distance between the primary source (PS) 204 and the exit surface of the object 202, and xa the distance between the exit surface and the x-ray receptor 206 (air gap). All the scattered x-ray photons could be assumed to originate from an effective scatter point source (ESPS) 208 at an xs distance from the exit surface of the object 202.
Governed by the inverse-square-law, both the primary and scattered x-ray flux are attenuated by an air gap. However, for a wide range of experimental conditions, such as the phantom thickness, FOV, x-ray source energy as well as the distance xp between the x-ray focal spot and the exit surface of the phantom, xs is consistently within the range of 15xcx9c20 cm and usually much shorter than xp. Consequently, the scattered x-ray photons are attenuated by an air gap to an extent larger than that of the primary x-ray photons. This means that, given an air gap within a reasonable range, more scattered x-ray photons are selectively attenuated by an air gap, thus improving SPR. A formula to get the optimal xa is                               x          a_opt                =                              1            2                    ⁢                      (                                          x                t                            -                              x                s                                      )                                              (        16        )            
where x1 is the allowable distance between the x-ray source focal spot and the receptor.
In conventional CT, the primary function of a bow-tie attenuator is to ameliorate the quality of the x-ray source, so that the artifact caused by beam hardening can be reduced. Another important function is to alleviate the requirement on the dynamic range of an x-ray detector and lower the radiation dose to a patient. Particularly, the SPR uniformity of a cylindrical object is improved, since the bow-tie attenuator compensates for the primary x-ray photon distribution non-uniformity, so that the cupping artifact in reconstructed images can be significantly reduced.
Another approach to suppress x-ray scatter is the anti-scatter grid, which is either linearly or crosswise focused. The linearly focused grid provides an inclusive performance better than that of the crosswise focused grid. In practice, the focal spot of an x-ray source is positioned on the convergent line of a linear grid that is perpendicular to the central beam emanating from the focal spot, so that most of the primary x-ray photons can transmit the interspaces. Since scattered x-ray photons behave as if they originated from its ESPS, most of them are attenuated by a linear grid. Generally, the larger the grid ratio, the more selective the attenuation. However, a grid with too high a ratio absorbs too many primary x-ray photons, resulting in an increased x-ray tube loading and exposure to the patient. Moreover, it was observed that the performance of both the air gap approach and the grid approach is strongly dependent on x-ray scatter. Whereas for high x-ray scatter both provide comparable improvement on SNR, for low to moderate situations an air gap definitely performs better. Furthermore, if x-ray scatter is very low(SPR less than 0.7), a grid even degrades SNR due to the absorption of primary x-ray photons by the grid. Therefore, for SPR less than 0.7, an antiscatter grid may not be used. For example, in cone beam volume CT breast imaging (CBVCT mammography) an antiscatter grid may not be used.
The other approach involves the use of a beam stop (BS) array. As shown in FIGS. 3A and 3B, an x-ray BS array 302 having small lead discs 304 is placed between the x-ray source 306 and the object to be imaged, which in an illustrative example is a phantom 308 on a table 310. On the far side of the phantom 308 from the source 306 is a detector 312. By adequately choosing the dimensions of the lead discs 304, each shadowed area in a projection image is as small as possible, while an entire block of the primary x-ray photons is assured. The small lead discs are deployed in an array with an orthogonal distance significantly larger than the diameter of small lead disc, so that the shadows are apart from each other. A total of two projection imagesxe2x80x94one with the beam stop array (image I) and the other without it (image II)xe2x80x94are taken. The intensity detected within each shadowed area in image II is assumed to be exclusively that of the scattered x-ray photons, resulting in the x-ray scatter distribution sampled on a sparse lattice. Employing a spatial cubic spline interpolation (or bilinear interpolation) on these samples, the scatter distribution (scatter image) is estimated. Subsequently, the scatter image is subtracted from image I to get the image exclusively formed by the primary x-ray photons (primary image).
The merit of the BS array approach is its ability of adaptively removing scattered x-ray intensity. Nevertheless, if it were directly employed in CBVCT, the resultant exposure would be doubled, and such a doubled exposure is definitely undesirable in practice.
The inventor""s previous work disclosed in WO 01/35829, published May 25, 2001, whose disclosure is hereby incorporated by reference in its entirety into the present disclosure, discloses a technique using a BS array in which the scaling factors for estimating scatter distribution and the convolution kernels at sampled angle positions can be determined. Then the scatter distributions are estimated using the convolution kernel at corresponding angle positions and subtracted from the detected positions. However, that technique requires determination of the convolution kernels, rather than the use of existing kernels for cubic spline interpolation. It also depends on compression of the patient""s breast into a cylindrical shape, which is not always desirable or even possible for any given object to be imaged.
It will be apparent from the above that a need exists in the art to improve the reconstruction accuracy in CBVCT without adversely affecting patient safety and without adversely reducing data acquisition speed. It is therefore an object of the present invention to reduce the x-ray scatter interference.
It is another object of the invention to reduce the x-ray scatter interference in CBVCT to a clinically acceptable level.
It is further object of the invention to reduce the x-ray scatter interference in CBVCT to improve reconstruction accuracy of CBVCT.
It is a still further object of the invention to reduce the x-ray scatter interference in CBVCT without the use of excessive additional x-ray exposure to the patient or, more generally, to the object to be imaged.
It is a still further object of the invention to reduce the number of images to be taken with the BS array without requiring the determination of convolution kernels at sampled angle positions.
It is a still further object of the invention to reduce the number of images to be taken with the BS array, regardless of the ability to compress the object to be imaged into any particular shape.
It is a still further object of the invention to combine the following techniques together to combat scatter interference in CBVCT: x-ray beam compensation filter (bow-tie filter), air gap technique, and antiscatter grid to reduce detected scatter, and an image processing technique for final scatter correction to improve reconstruction accuracy.
To achieve the above and other objects, the present invention is directed to a system and method for imaging in which the x-ray scatter interference is reduced through the decomposition of primary and scatter image subsequences. From the perspective of image sequence processing, all the consecutively acquired projection images constitute an image sequence that includes a primary image subsequence (PIS) needed by CB reconstruction and a scatter image subsequence (SIS) interfering with CB reconstruction. The inventor has found that the spatial fluctuation of the SIS is significantly less intense than that of the PIS, indicating that a considerable spatial correlation exists in the former. Hence, it is possible to recover an SIS from samples acquired at a very low sampling rate using the BS array technique, and then its associated PIS can be obtained by correspondingly subtracting the SIS from the projection image sequence. The sampling rate of an SIS could be so low, e.g., 1, to 2/rotation 4/rotation, that its samples can be obtained in a way similar to the acquisition of scout localization images, which is usually employed in conventional/spiral CT to prescribe scan location relative to anatomic landmarks. Consequently, the resultant extra exposure and extra acquisition time would be acceptable. Again, the cubic spline interpolation or bilinear interpolation is angularly used to recover an SIS from its samples.
Theoretical and experimental investigations conducted by the inventor reveal that the anti-scatter approaches inherited from projection imaging, such as the air gap, beam-shaping attenuator and linear anti-scatter grid, can substantially reduce cupping/shading distortion (LAC reconstruction inaccuracy) in a CBVCT image caused by x-ray scatter interference, and so does the present invention. The quantitative performance evaluation of the present invention carried out using the scatter phantom shows that its capability of adaptively suppressing x-ray scatter interference in CBVCT imaging is still satisfactory even at an AI (angular interval) of 90xc2x0, which is equivalent to the acquisition of only 4 extra projection images with the BS array installed. This means that the present invention provides a novel and effective approach to combat x-ray scatter interference in CBVCT, and the extra x-ray exposure induced by it is definitely tolerable in practice.
A particular embodiment of the present invention provides a hybrid method by optimally combining the following techniques (discussed above) together to combat scatter interference: an x-ray beam compensation filter (bow-tie filter), an air gap technique, an antiscatter grid to reduce detected scatter, and an image processing technique for final scatter correction to improve reconstruction accuracy.
Finally, the penumbra phenomena resulting from the finite focal spot of an x-ray tube disturbs an accurate sampling of the distribution formed by scattered x-ray photons. Calibration can be used in practice to alleviate, if not eliminate, such a disturbance as much as possible.